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ABSTRACT 

The power spectram of density fluctuations is a foundational source of cosmological information. Precision 
cosmological probes targeted primarily at investigations of dark energy require accurate theoretical determinations 
of the power spectrum in the nonlinear regime. To exploit the observational power of future cosmological surveys, 
accuracy demands on the theory are at the one percent level or better Numerical simulations are currently the 
only way to produce sufficiently error-controlled predictions for the power spectrum. The very high computational 
cost of (precision) A^-body simulations is a major obstacle to obtaining predictions in the nonlinear regime, while 
scanning over cosmological parameters. Near-future observations, however, are likely to provide a meaningful 
constraint only on constant dark energy equation of state ' wCDM' cosmologies. In this paper we demonstrate that 
a limited set of only 37 cosmological models - the "Coyote Universe" suite - can be used to predict the nonlinear 
matter power spectrum at the required accuracy over a prior parameter range set by cosmic microwave background 
observations. This paper is the second in a series of three, with the final aim to provide a high-accuracy prediction 
scheme for the nonlinear matter power spectrum for wCDM cosmologies. 

Subject headings: methods: statistical — cosmology: large-scale structure of the universe 



1. INTRODUCTION 



Although the disco very of cosmic acceleration bv lRiess et al.l 
(11998) and Perlmutte r et al.l ( 119991) is already a decade in the 
past, our understanding of the nature of the underlying driver of 
the acceleration, "dark energy", has made little progress. One 
reason for this is that the dark energy equation of state param- 
eter w is consistent with a cosmological constant (w = -1) at 
roughly 10% accuracy, with no constraints on any possible time 
dependence. In order to advance further in terms of distinguish- 
ing different models of dark energy from each other and dark 
energy itself from other possible causes of acceleration, such as 
a possible break-down of general relativity on very large scales, 
observational errors must be brought down significantly. The 
current target is to achieve another order of magnitude improve- 
ment for several dark energy probes - probes that measure not 
only the expansion history of the Universe but also the growth 
of cosmological structure - down to the level of a percent. 

To date, the four most promising lines of investigation are: 
(i) Supernovae Type la, to measure the expansion history of the 
Universe, (ii) clusters of galaxies, to measure the expansion his- 
tory and growth of structure, (iii) baryon acoustic oscillations, 
to measure the expansion history, and (iv) weak lensing, to mea- 
sure the expansion history and the growth of structure. The last 
two probes, baryon acoustic oscillations and weak lensing, rely 
the most strongly on precise predictions of the nonlinear matter 
power spectrum. Numerical simulations are essential for cary- 
ing out this task, not only for the power spectrum itself, but 
also to build the underlying skeleton of cosmological structure 
from which object catalogs can be constructed. The resulting 
'mock catalogs' have many uses: to design and test observa- 
tional strategies, to understand systematic errors therein, and to 
confront theoretical predictions with observations. 



In the case of baryon acoustic oscillations, measurements are 
carried out on very large scales, where the nonlinear effects are 
small. Therefore, higher order perturbation theory might offer 
an alternative path to obtaining precise predictions for the non- 
linear matter power spectrum (see, e.g .,jCrocce & Scoccim arrol 
2006; Matsubara 2008; Pietronill20oi ICarlson et al.l i2009 and 
references therein), and provide a useful foil for the numerical 
results. Weak lensing measurements go down to much smaller 
spatial scales, out to wavenumbers k^^ 1-10 /zMpc"' (and even 
larger wavenumbers in the future). On these smaller length 
scales, perturbative techniques fail, and one must rely on nu- 
merical simulations to obtain the required level of accuracy: at 
least as accurate as the observations, and t o be optimal, sub- 
stantially more accurate. As shown by, e.g. iHuterer & Takadal 
(2003), to avoid biasing of cosmological parameter estimations 
a wide-field weak lensing survey such as the SNAld design re- 
quires 1 % accurate power spectrum predictions, and a survey 
such as the Large Synoptic Survey Telescope (LSSlQ) requires 
predictions at the 0.5% accuracy level. 

These requirements pose two major challenges: First, one 
must show that simulations capturing the essential physics have 
reached the desired level of accuracy. For baryon acoustic os- 
cillations, it is expected that gravity-only A^-body simulations, 
augmented by halo occupancy modeling, are sufficient for the 
task. In the case of weak lensing, this assumption holds for 
scal es out to ^ ^ 1 /zMp c"'. In the first paper of this se- 
ries (iHeitmann et al.l2008l) we have established that, up to these 
scales, the nonlinear matter power spectrum can be determined 
at sub-percent accuracy by gravity-only A^-body simulations. 
At smaller scales, baryonic physics becomes important at the 
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few to ten percent level and has to be taken int o accoun t ( White! 
l2004tlZhan & Knoxl2004l:rjing et al. 2006; Ru dd et al.ll2008h . a 
task which has to be tackled accurately in the near future, per- 
haps by a suitable combination of simulations and observations. 

After overcoming the first challenge, the next task in con- 
straining cosmological parameters, is to cover a range of differ- 
ent cosmologies. Markov Chain Monte Carlo (MCMC) meth- 
ods, commonly used for parameter determination, rely on re- 
sults from model evaluations numbering in the tens of thou- 
sands to hundreds of thousands. Since an accurate A^-body sim- 
ulation on the scales of interest mentioned above costs of the 
order of ~ 20,000 CPU-hours, it is not feasible to run such 
simulations for each model. (Running ^ 20,000 A^-body sim- 
ulations with the required resolution on a contemporary 2048 
processor cluster would take 20 years!) Taking into account the 
fact that adding gasdynamics and feedback effects substantially 
increases both the computational load and the number of pa- 
rameters to be varied, it is clear that a brute force approach to 
the problem has to be avoided. What we need is a generalized 
interpolation method capable of yielding very accurate predic- 
tions for the nonlinear matter power spectrum from a restricted 
number of simulations. In the following, we will refer to such a 
prediction scheme as an emulator The emulator will be tasked 
with replacing brute force A^-body simulations for the nonlinear 
matter power spectrum over a pre-defined set of cosmological 
parameters, with specified ranges for the chosen parameters. 

In the cosmic microwave background (CMB) community 
several different paths have been suggested to provide such an 
emulator for the CMB temperature anisotropy power spectrum. 
These include purely analytic fits (Tegmark & Zaldarriaga 2000; 
iJimenez et al . 2004) and cornbinat ions of analytic and semi- 
analytic fits tKaplinghat et al.l2002l) . More recently, neural net- 
work methods and machine learning techniques have been suc- 
cessfully used to generate very accurate temperature anisotropy 
power spectra (Fendt & Wandelt 2007a; Auld et al. 2007; Fendt 
& Wandelt 2007b; Auld et al. 2008). These methods are based 
on a large number of training sets, up to several tens of thou- 
sands. While this does not constitute a problem for anisotropy 
power spectra - given the speeds at which numerical codes such 
as CAME and CMBFast can be run - the approach is not feasi- 
ble for determining the matter power spectrum, which requires 
large-scale supercomputer simulations. 

As in the case for the temperature anisotropy power spec- 
trum, several attempts have been made to avoid costly simula- 
tions by finding good approximations for the nonlinear matter 
power spect rum. These range from more or less analytic deri va- 
tions (e.g., [Hamilton etalTi 19911; iPeacock & Doddsl Il994l) to 
semi-analytic fits calibrated more specifically against simula - 
tion results (e.g.. Peacoc k & Doddslll996l; ISmith et al]|2003l) . 
Unfortunately, the accuracy of these approxim ations is inade- 
quate , at best reaching the 5-10% level (see, e.g.. lHeitmann et all 
|2008| for a recent comparison of precision s imulations with 
HaloFit, a fitting scheme due to Smith et a l. 2003). Thus, 
an order of magnitude improvement is needed to address the 
accuracy requirements discussed above. 

Accurate emulation is needed for many observational quan- 
tities in cosmology, power spectra being one important ex- 
ample. To address this need, we ha ve recently introduced 
the "Cosmic Calibration Framework" (jHeitmann et aT] l2006t 
iHabibet all 120071: [Schneider et all l2008l) combining sophisti- 
cated simulation designs with Gaussian process (GP) model- 
ing to create very accurate emulators from a restricted set of 



simulations. The term 'simulation design' refers to the spe- 
cific choice of parameter settings at which to carry out the 
simulations. One of the main reasons why the Cosmic Cali- 
bration Framework provides very accurate results from only a 
small number of training sets is the optimization of the simula- 
tion design strategy to work well with the chosen interpolation 
scheme, the Gaussian process in this particular case. Another 
useful aspect of the methodology is that it contains an error 
prediction scheme, so that one can verify the consistency of the 
obtained results. 

In this paper we will explain and demonstrate the emulation 
capability of the Cosmic Calibration Framework. With only a 
small number of simulations, an emulator for the nonlinear mat- 
ter power spectrum can be constructed, matching the simulation 
results at the level of 1 % accuracy. We focus on the regime of 
spatial wavenumber k<l /iMpc"' and a redshift range between 
z = and z = 1, covering the current space of interest for weak 
lensing measurements. Such an emulator will eliminate a ma- 
jor source of bias in determining cosmological parameters from 
weak lensing data. In order to design, construct, and test an em- 
ulator, it is useful to carry out the process first on a proxy for the 
expensive numerical simulations; the proxy need not be very 
accurate but should reflect the overall behavior of the detailed 
simulations; we employ HaloFit in this role. 

This paper is the second in a series of three communications. 
In the first, we have demonstrated that it is possible to obtain 
nonlinear matter power spectra at sub-percent level accuracy 
out to A; ^ 1 /iMpc"' from simulations, having derived and pre- 
sented a set of requirements for such simulations. The third 
paper of the series will present results from the complete sim- 
ulation suite based on the cosmologies presented in the current 
paper, as well as the public release of a precision power spec- 
trum emulator. The simulation suite is named the "Coyote Uni- 
verse" after the computer cluster it has been carried out on. 

The paper is organized as follows. In Section |2] we describe 
in detail the Cosmic Calibration Framework with special em- 
phasis on building a nonlinear matter power spectrum emulator 
from a restricted set of simulations. We explain the design strat- 
egy for generating the training sets and discuss the emulation 
step, demonstrating the emulator accuracy. Next we provide a 
sensitivity analysis showing how the power spectrum varies - 
in a scale-dependent manner - as the cosmological parameters 
are changed. We conclude in Section[3] 

2. THE COSMIC CALIBRATION FRAMEWORK 



The Cosmic Calibration Framework (iHeitmann et al.l 120061; 
iHabib et al.ll2007HSchneider et al.ll2068l) consists of four inter- 
locking steps: (i) the simulation design, which determines at 
what parameter settings to generate the training sets, (ii) gener- 
ation of the emulator which replaces the simulator as a predictor 
of results away from the points that were used to generate the 
training set, (iii) the uncertainty and sensitivity analysis associ- 
ated with the emulator, and (iv) the calibration against data via 
MCMC methods to determine parameter constraints. 

In the following we discuss steps (i) - (iii) in detail, with spe- 
cial emphasis on generating an accurate emulator for the non- 
linear matter power spectrum. 

2.1. Sampling the Model Space 

As discussed in the Introduction, one of the major challenges 
in building an accurate emulator for the nonlinear matter power 
spectrum is the very high cost of individual A^-body simula- 
tions combined with the high dimensionality of the parameter 



Heitmann, Higdon, White, Habib, Williams, Wagner 



3 



space (which may include cosmological, physical, and numer- 
ical modeling parameters). To sample the model space, the 
number of parameters to be varied must be specified, as well 
as the range of variation for each parameter For now, we will 
assume that some combination of observational knowledge and 
cosmological and astrophysical modeling is sufficient to decide 
on conservative choices for sampling model space. (We will 
return to this question later, in Section l2.1.2l ) 

Following this decision, the next step is to find a method for 
sampling the model space and interpolating the results there- 
from, satisfying given accuracy criteria, and using only a man- 
ageable number of simulation design points. In several appli- 
cations, space-filling Latin hypercube (LH) designs (details be- 
low) ha ve proven to be well suited for the GP model-based ap- 
proach dSacks et alJl989atlCurrin et al.lll991h to solving the in- 
terpolation problem. The Cosmic Calibration Framework uses 
this met hodology ; the associated validation examples can be 
found in lHeitmann et al] (l2006l) and lHabib et al.1 UQQlh . 

In the following, we first discuss the statistical aspects of 
sampling model space followed by the observational aspects, 
i.e., the parameters to be varied and the corresponding choices 
for the range of variation. The parameter choices and prior 
ranges from observations used he re rely on the most recent 
CMB observations from WMAP-5 dKomatsu et al.ll2008l) . 

2.1.1. Statistical Sampling Methods 

Our first aim is to find a distribution of the parameter settings 
- the simulation design - which provides optimal coverage (in 
a sense to be defined below) of the parameter space, using only 
a limited number of sampling points. (In the statistics literature, 
it is customary to use normalized units in which all parameters 
range over the interval [0, 1] and we will follow this usage for 
the most part.) If the actual behavior of the observable as a 
function of the parameters is considered to be unknown, then 
it is sensible to start with a strategy that attempts to uniformly 
sample the space (space-filling design). An extreme version 
of this is a simple, regular hypercubic grid. The problem with 
using a grid-based interpolation method is the high dimension- 
ality of the space. Suppose we wish to vary 5 cosmological pa- 
rameters and sample each parameter only three times - at near- 
maximum, near-minimum, and in the middle of each parameter 
range. Already, this would require 3^ = 243 simulations - not 
a small number - and lead to poor coverage of the parameter 
space. Such a design with only 3 levels (three sampling points 
per parameter) would also only allow for estimating a quadratic 
model. If we want to go to a higher level, the number of runs 
will explode - if we wish to keep the complete grid. If on the 
other hand, we try to reduce the number of runs by using only 
a fraction of the grid, aliasing becomes a potential problem. 

The opposite extreme of pure random sampling suffers from 
clustering of sample points and occurence of large voids in the 
sampling region when the number of sample points is limited. 
Stratifed sampling techniques combine the idea of a regular grid 
and random sampling by using strata that (equally) sub-divide 
the sample space, with random sampling employed within each 
strata. A final important point is that the computed observable 
may depend on some sub-combination of input variables (factor 
sparsity), in which case we would like to have uniform coverage 
across the projection of the full space onto the relevant factors. 
Not all uniform sampling schemes possess this property. 

The GP model interpolation scheme used here does not re- 
quire a simple grid design. Simulation designs well-suited for 



GP model emulators are LH-based designs, a type of stratified 
sampling scheme. Latin hypercube sampling generalizes the 
Latin square for two variables, where only one sampling point 
can exist in each row and each column. A Latin hypercube sam- 
ple - in arbitrary dimensions - consists of points stratified in 
each (axis-oriented) projection. More formally, an LH design 
is an n X m matrix in which each column is a unique random 
permutation of {1, ...,n}. The use of LH designs in applications 
where the aim is to predict the outcome of some quantity at un- 
tried parameter se ttings from a restrict ed set of simulations was 
first proposed by ' McKay et al.l (II979'). As discussed in more 
detail below, like many other interpolators, GP models rely on 
local information for their interpolation strategy. Therefore, the 
simulation design must provide good coverage over the whole 
parameter space. Space-filling LH designs and their variants 
provide an optimal approach for achieving this. 

Very often LH designs are combined with other design strate- 
gies such as orthogonal array (OA)-based designs or are opti- 
mized in other ways, e.g., by symmetrizing them (more details 
below). By intelligently melding design strategies, different at- 
tributes of the individual sampling strategies can be combined 
to lead to improved designs, and shortcomings of specific de- 
signs can be eliminated. As a last step, optimization schemes 
are often applied to spread out the points evenly in a projected 
space. One such optimization scheme is based on minimizing 
the maximal distance between points in the parameter space, 
which will lead to more even coverage. For a discussion of dif- 
ferent design approa ches and their specific advantages see, e.g., 
ISantner et al](l2003h . 

We now briefly discuss two design strategies well suited to 
cosmological applications in which the number of parameters 
is much less than the number of simulations that can be per- 
formed. These are optimal OA-LH design strategies and op- 
timal symmetric LH design strategies. The former has been 
used in previous work in cosmology (Heitmann et al. 2006; 
iHabib et al. while the latter will be used in this paper 

to construct the design for the Coyote Universe. For illustration 
purposes, we will use a very simple, three-dimensional case 
example with three parameters, 6*1, 02, O3 and nine sampling 
points. 

In order to create an OA-LH design, the strategy proceeds 
in two steps: (i) construction of the orthogonal array and (ii) 
the following construction of the orthogonal-array based Latin 
hypercube. We discuss these two steps in turn following closely 
the description bv , Tang. C 1993.) and .Learv et al.. (.20031) . 

OA Designs 

An orthogonal array distributes runs uniformly in certain pro- 
jections of the full parameter space. The mathematical defini- 
tion is as follows: An n by m matrix A with entries from a set 
1 , 2 , .... 5 is called an orthogonal array of strength r, size n, with 
m constraints, and s levels, if each nx r submatrix of A con- 
tains all possible 1 x r rows with the same frequency A. Here 
A is termed the index of the array , and n = Xs''. The array is 
denoted OA(n, m, s, r) (lTanglll993h . 

For our application, n denotes the number of simulation runs 
to be performed and m specifies the number of parameters to be 
varied (these can be cosmological parameters as well as numer- 
ical input parameters). These choices fix the number of dimen- 
sions in the parameter hypercube. The parameter s defines the 
levels of stratification for each column in the matrix A. In order 
to sample the parameter space well in a uniform manner, it is of- 
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ten not enough to sample it well globally. For example, if two or 
more parameters interact strongly with each other, it is clearly 
desirable to have a good space-filling design in the subspace 
of these parameters. In other words, if one projects the design 
down onto, e.g, two dimensions, such a projection should have 
space-filling properties in those dimensions as well. The pa- 
rameter r, the strength of OA designs, indicates the projections 
for which the LH design based on that OA are guaranteed to be 
space-filling. For example, if r = 3, then all 1,2 and 3 dimen- 
sional projections will be space-filling. Obviously, r cannot be 
larger than m, the number of parameters varied. 

The strength r is usually restricted to two or three for several 
reasons: (i) Fewer and fewer known OA designs with appropri- 
ate run sizes exist as the strength increases, (ii) In most applica- 
tions, only a small number of parameters influence the response 
significantly. Statisticians call this the "20-80 rule" - 20% of 
the parameters being responsible for 80% of the outcome varia- 
tion. Therefore, the aim is to capture the action of these relevant 
parameters. Furthermore, outcome variation is often dominated 
by a small number of single parameter and two-factor interac- 
tion effects, which are adequately covered by OA-LH designs 
based on r = 2 or 3. (iii) The number of simulations often has to 
be kept small, therefore r cannot be chosen too large, since the 
number of simulations n is connected to r via n = \s'' . As is the 
case for r, the stratification parameter s is also restricted by the 
number of runs one can possibly perform. It is very often set to 
s = 2 which then requires the number of runs to be a power of 
two. The frequency A with which a permutation repeats, is kept 
small as well. To create a design, the strategy is to fix strength 
first, and try to find an OA design that has approximately the 
right number of runs and at least as many parameters as one 
needs. If such a design cannot be found, then the strength is 
reduced by one and the process repeated. Usually, this strategy 
is started with OA designs of strength 3, though many more de- 
signs of strength 2 exist. It is rarely possible to find a strength 
four or higher design with few enough runs. 

The above discussion points to a shortcoming of orthogonal 
arrays: the number of simulation runs cannot be picked arbi- 
trarily (e.g., choosing s = 2 forces a power of two for the num- 
ber of runs). The other difficulty with orthogonal arrays is that 
they are not easy to construct. A detailed desc ription of OA 
design s and how to construct them is given in iHedavat et al.l 
(Il999l) . A library containing a large number of OAs can be 
found at: http://www.research.att. com/^njas/oadir/index.html 
A collection of C routines to generate OA designs can be found 
at: http://lib.stat.cmu. edu/designs/oa.c 

We now turn to our example with nine sampling points (« = 9) 




Fig. 1. — Left panel: an orthogonal array (OA) based design for 3 param- 
eters, 9i, 02, 03 and nine sampling points. Right panel: the OA based design 
perturbed in such a way that the one-dimensional projection onto any param- 
eter leads to an equally spaced distribution of sample points. The projection 
onto any two dimensions leads to a space filling design. 
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Fig. 2. — Projections of the design shown in Figure[T]onto two dimensions. 
The lower triangle shows the projection of the OA design, the upper triangle of 
the OA-LH design. Note that when projected onto one dimension, the OA-LH 
design leads to an even coverage and no points lie on top of each other. 



and three parameters (m = 3). If we require r > I for the 
strength of the design, then the relation n = Xs'' leads automati- 
cally to an OA design with s = 3 levels, strength r = 2, A = 1 ; an 
OA(9,3,3,2) design. We require that if we project our design 
down onto any two-dimensional direction, the parameter space 
be well covered. The left panel in Figure [T] shows a possible re- 
alization of an OA with our example specifications. The lower 
triangle in Figure |2] shows the three possible two-dimensional 
projections of this design. This specific design is of course not 
a unique solution. In matrix form it reads: 
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(1) 



From this 9x3 matrix we can now verify that each of the 
three 9x2 sub-matrices indeed contain all possible 1x2 rows 
with the same frequency A = 1 . On the right hand side of Eq. ([T]) 
we have simply rescaled the design points into the normalized 
[0,1] space which is shown in Figure s[T| and |2] 

OA-LH Designs 

In order to further improve parameter space coverage, the next 
step - Latin hypercube sampling - perturbs the positions of 
each sampling point from A via the following prescription: for 
each column of A, the Ai*^"' positions with entry k are replaced 
by a permutation of (k= 1 , • • • ,s) 



Heitmann, Higdon, White, Habib, Williams, Wagner 



5 



(k- 1)A.?'-' + l,(k- l)X/-^ + 2,--- ,(k- l)As'-' + X/-^ = kX/-^ 

(2) 

This means, in our example, that the entries for k = I will 
be replaced by 1,2,3, the entries for k = 2 will be replaced by 
4,5,6, and the entries for A: = 3 by 7,8,9. The Latin hypercube 
algorithm demands that in every column every entry appears 
only once. This ensures that each one dimensional projection 
is evenly covered with points and no run is replicated in the 
resulting design. The right panel in Figure [T] shows a possible 
realization of this in three dimensions, derived from perturbing 
the orthogonal array in the left panel. The upper right triangle 
in Figure |2] shows the two-dimensional projection of this de- 
sign. The solid lines show the division for the orthogonal array 
while the dashed lines show the additional sub-division. Note 
that each sample point lies on a unique horizontal and verti- 
cal dashed line - if the sample points are projected down onto 
any one direction, the one-dimensional space would be evenly 
covered. In matrix form, our OA-LH design is as follows: 
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Note that we have replaced the entries in a random fashion 
and created this design "by hand" - convincing o urselves "by 
eye" th at we have good coverage in 2-D projection. iLearv et al] 
( l2003h suggest an optimal strategy to ensure even better cover- 
age of the parameter space. These optimization strategies are 
often used for the projected parameter space. In order for the 
points to be spread out, one has to determine the 'closeness' 
between them. This can be defined as the smallest distance be- 
tween any two points. A design that maximizes this measure 
is said to be a max imin distance design. ( For more details , see 
Santner et al. II2003L ') The designs m Heit mann et al. I (l2006l) and 



Habib et aU (l2007h combine the OA-LH based design with a 



maximin distance design in each two-dimensional projection. 
Other optimization methods rest on an entropy criterion based 
on the minimizati on of -log \R\, where R i s the covariance ma- 



trix of the design (IShewrv & Wvnnl 1987h. or rn inimization of 



the Integrated Mean Squared Error (I Sacks et al ] [T989bl) . 

Our example shows just one way to realize an OA-LH de- 
sign. It can be implemented straightforwardly and leads to the 
desired coverage of the parameter space. After the design has 
been determined in the [0,1] parameter space, it then can be eas- 
ily translated into the physical parameter space of interest. At 
this point, when projected down to one dimension, the equidis- 
tant coverage in each dimension of the parameter space in one 
dimension is of course lost. However, since our transforma- 
tion is linear, we do not lose the uniformity of the projections. 
Therefore it is still true that no two sample points will fall on 
top of each other in projection. 

SLH designs 

As mentioned above, the major restriction of OA-LH based de- 
signs is that they cannot be set up for an arbitrary number of 



simulation runs. This is a specific point of concern, if one can 
only run a very restricted number of simulations. In addition, 
the set-up of an OA-LH design can be nontrivial. Very often, 
one has to rely on OA libraries which are restricted in their pa- 
rameters and also not always easily availa ble. An altern ative 
space-filling design strategy presented by Li & Yd (l2000l) . of- 
fers a compromise between the computing effort and the design 
optimality - the optimal symmetric Latin hypercube designs. 
Following their definition, an LH design is called a symmetric 
LH (SLH) design if it has the following property: For any row 
/ of an LH design, there exists another row in the design which 
is the /th row's reflection through the center. For example, in an 
nxm Latin hypercube with levels from 1 to n, if (01,02, ■■■,am) 
is one of the rows, the vector (n+ 1 -oi,n+l -02, ...,n+ l-Om) 
should be another row in the design. The symmetry imposes 
a space-filling requirement on the designs considered up front, 
which carries through to all projections. An example for an 
SLH design is given by: 



(4) 



In this design, rows 1/2, 3/4, 5/6, and 7/8 are symmetric pairs. 
As before, we do not attempt to optimize the resulting design, 
though the design we will use for the remainder of the paper 
will be optimzed, see below. The two-dimensional projection 
of the design given in Eq. (|4|i is shown in Figure[3] 

iLi & Ye (2000) provide an excellent discussion of optimal 
SLH designs, including a description of possible algorithmic 
implementations and comparison with traditional optimal LH 
designs. As an example, they show that the computational ef- 
fort to find an optimal LH design by starting with an SLH de- 
sign reduces roughly by a factor of ten for a 25 x 4 design on 
a workstation. As before, the SLH design is usually optimized 
in the last step, often with respect to a distance based criterion 
which spreads out the points in two-dimensional projections. 
Two standard search algorithms for an optimal S LH de sign 
are the columnwise-pairwise (CP) algorith m bv,"^ ( 19981) and 
the sim ulated annealing (S A) algorithm by iMorris & Mitchelj 
( 1 19951) . Simply put, these algorithms are based on columnwise 
exchanges of entries which will keep the symmetry properties 
(since the corresponding symmetric pairs are always switched 
at the same time). They are iterative procedures, which will 
stop after a certain pre-set optimization criterion is fulfilled or 
the process is interrupted by time limitations. Very often, sev- 
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Fig. 3. — Two-dimensional projections of tlie SLH design given in Eq. j4). 
The symmetric design points are connected to sliow tlie reflection through the 
center. 
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eral designs are produced at the same time and the most optimal 
kept in the end. The SA and CP algorithms can also be used to 
optimize OA-LH designs. If the OA skeleton is symmetric, then 
one can require the optimal OA-LH design to be symmetric as 
well. 

In the following, we will use an LH design optimized via a 
distance criterion. The design will encompass 37 simulation 
runs and five cosmological parameters. In detail, 20 optimiza- 
tions with CP and 20 with S A were carried out, and the best was 
chosen in the end where the quality was measured by a distance 
criterion. For CP, 10 of the designs had a symmetry require- 
ment and the other 10 did not. For SA, 10 of the designs had 
a symmetry requirement and were optimized with a local op- 
timization criterion, and the other 10 did not have a symmetry 
requirement and were optimized with a more global optimiza- 
tion criterion. The best design from all of these came from one 
of the optimizations using SA, a global optimization criterion, 
and no symmetry requirement. 

2.1.2. Observational Considerations 

The choice of number of active parameters depends on the 
available data as well as on the chosen modeling approach. We 
do not insist on a formal methodology to make this choice, but 
instead present practical and conservative arguments to justify 
our decisions. We take as our basic 5 parameters uj,„ = Vl,„h^, 
ujb = ^bh^, ns, w, and erg where il,„ contains the contributions 
from the dark matter and the baryons. We restrict ourselves to 
power-law models (no running of the spectral index), to spa- 
tially flat models without massive neutrinos and to dark energy 
models with constant equation of state. A sixth parameter, the 




400 600 
Multipole Moment / 



1000 



Fig. 4. — Best-fit TT power spectra for each model in Table [T] using the 
WMAP-5 results. The only parameter which has been optimized by mini- 
mizing is the optical depth r. The upper panel shows the resulting power 
spectra, the black points with error bars show WMAP-5 data points, and the 
thick black line the best-fit WMAP-5 model. The lower panel shows the resid- 
uals for each model with respect to the best-fit model. Some of our models are 
undemormalized, the best-fit t being smaller than 0.01 which would lead to a 
reionization redshift of z < 2 and x for these models is lai'ger than 3000 (the 

for the best-fit WMAP-5 model is at roughly 265' 
models at r = 0.01 and show them with dashed lines. 



redshift or time, simply requires us to dump data from each run 
at multiple epochs. 

The effect of massive neutrinos can be roughly approximated 
by decreasing fl,,, (Brandbyge et al. 2008). We do not ex- 
pect any realistic dark energy model to have a constant equa- 
tion of state, but we wanted to begin with the most restric- 
tive parameter space in order to validate our methods. The 
next generation of experiments will pose at best weak con- 
straints on any time variation of w, and in this sense our 
constant w can be thought of as an ap propriate average of 
w(z). Usin g growth matching techniques (I White & Vald 120041 : 
iLinder & W hite 2005; Francis et al. 2007) one can map the 
power spectrum from a complex w{z) onto one with a constant 
w with reasonable accuracy. 

The normalization is another area where choices need to be 
made. Historically the amplitude of the power spectrum was set 
by (78, the amplitude of the linear theory matter power spectrum 
smoothed with a top-hat on scales of 8 /!"'Mpc 



dk . T 



3jdkR) 



kR 



«=8/7-iMpc 

with the linear power spectrum being defined as 

k^PiJk) 



27r2 



(5) 



(6) 



This scale and normalization were chosen because the fluctu- 
ations of counts of galaxies in cells of this size is close to 
unity. With the advent of the COBE data it bec ame common to 
quote the normalization at horizon scales, e.g. iBunn & Whitj 
( Il997h . As CMB data improved, the pivot point was shifted 
to smaller scales, closer to the middle of the range over which 
the spectrum is probed and where the normalization is best de- 
termined. In order to make closer connection with the initial 
fluctuations, the amplitude not of the matter power spectrum 
but of the curvature or potential fluctuations has been adopted. 
These differ mostly by factors of growth and fl,„. Anticipat- 
ing future advances kp = 0.002Mpc~' was sele cted for the most 
recent CMB analysis bv lKomatsu et al.l (l2008l) and the rms cur- 
vature fluctuation on this scale is now most commonly used 
as a normalization. With present CMB data the biggest uncer- 
tainties in the normalization are the near degeneracy with the 
optical depth, r, and the uncertain growth of perturbations at 
low redshif t due to the u nknown equation of state of the dark 
energy, e.g. lWhitel(l2006l) . 

For our purposes, however, a normalization tied to the present 
day matter power spectrum and close to the nonlinear scale has 
many advantages. Rather than introduce yet another conven- 
tion, we therefore choose to use erg as our normalization pa- 
rameter. Of course, since all of the parameters of the models 
are specified one can compute any other parameter for our mod- 
els. As an example, we have evaluated for each of the models 
given in Table [1] the best-fit value for r using the likelihood 
code provided by the WMAP-5 team. The resulting TT power 
spectra are shown in Figure |4] as well as their ratios with re- 
spect to the best-fit WMAP-5 model. Some of our models are 
undemormalized and the resulting r is smaller than 0.01 lead- 
ing to reionization redshifts of z < 2. This undernormalization 
however is not a concern: we chose the 38 models to cover the 
parameter space well overall and not to provide fits close to the 
concordance cosmology. We provide the best-fit values for r in 
Table ID 

From the WMAP 5 -year data, in combination with BAO, we 
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hav^l 



tj,„= 0.1 347 ±0.0040 (3%), 
a;fo = 0.0227 ±0.0006 (3%), 
n, = 0.9610 ±0.0140 (2%). 



(7) 



Current data restrict a constant equation of state for the dark en- 
ergy to w = -1 with roughly 10% accura cy (for ver y recent re- 
sults from supern ovae see, e.g.,,Kowalski et al. I I2OO8: for weak 
lensing see, e.g., lKilbingeretal.1120081; and for the latest con- 
straints from clusters of galaxies, see lVikhlinin et al.l2008l) . Re- 
cent determinations put the normalization in the range 0.7-0.9 
with still rather large uncertainti es (see, e.g. , Vikhlinin et al. 
l2008l for constraints from clusters, Voevodkin & yikhlinia2004i 
for a low estimate from clusters, iTegmark et al.l 2007 for con- 
straints from combined CMB and large scale structure data, 
and Evrard et al. 2008 for an extended discussion of recent re- 
sults). Considering all these constraints and their uncertainties, 
we choose our sample space boundaries to be 



0.120 <tJ,„< 0.155 
0.0215 <cjfc< 0.0235 
0.85 <«5< 1-05 
-1.30< w<-0.70 
0.61 < (78 < 0.9 



(8) 



as shown in Table [T] 

In order to solve for the full set of cosmological parame- 
ters we impose the CMB constraint that = Trt/is/r? = 302.4, 
where c/is is the distance to the last scattering surface and is 
the sound horizon. Observationally this is known to 0.3%, and 
models which significantly violate this equality are poor fits to 
the CMB data (see Figure|5]). Unfortunately the sound horizon, 
like the epoch of last scattering, can be defined in a number of 
different ways which differ subtly. Specific ally we use the fit to 
the redshift of last s cattering of Eg. (23) ofjHu & White! ( Il997h 
and use Eq. (B6) of Hu & Sugivama ( 1995) for the sound hori- 
zon. With these choices we find models with lli,„ and ujh in the 
range preferred by WMAP do indeed provide good fits to the 
WMAP data. This is demonstrated for model 32 in Figure|5] 

The procedure is then as follows. For every specified lj,„ 
and LOb we compute and hence the required to fit the £a 
constraint. We adjust h, at fixed spatial curvature, w, and ujm, 

^ See [http://lambda.gsfc.nasa.gov/ 1 




0.75 0.8 0.85 

Hubble parameter h 

Fig. 5. — Sweep thi'ough h for model 32. The red circle marks the estimate 
for the Hubble parameter from assuming perfect knowledge of Ia , in excellent 
agreement with the result from the WMAP-5 likelihood for the best-fit value 
of h for this model. 



until the model reproduces the required t/k. Knowing h and uj,„ 
then gives us il,„ and hence fide, as shown in Table |2] 

Finally, we generated a model '0' which has parameters 
close to t he current best fit fro m CMB and large-scale struc- 
ture (e.g.. iKomatsu et aDl2008h . This model has f7„, = 0.25, 
f^de = 0.75, ujb = 0.0224, n, = 0.97, h = 0.72, and = 0.8 and 
can be used as an independent check of the interpolation in the 
range of most interest. 

2.1.3. The Resulting Design 

Based on the above considerations, we can now generate a 
space-filling design for the five parameters of interest. We re- 
strict ourselves to 37 cosmologies and will show in the remain- 
der of the paper that this number is indeed sufficient to generate 
an accurate emulator. For model we pick a standard LCDM 
model for which we chose the Hubble parameter by hand (al- 
though h = 0.72 is very close to the result we would obtain if 
we would derive it from t/^). The resulting cosmological mod- 
els are listed in Table [T] where we give the values for the basic 
parameters. In Table |2] we give in addition a few derived pa- 
rameters: fim, f^de (recall that flatness is assumed), h as derived 
from our constraint on £a, d\s, and r. 

The two-dimensional projection of the design is shown in 
Figure |6] The upper part of the triangle shows the five input 
parameters in red, demonstrating a good sampling of the pa- 
rameter space. The blue points show projections of three of the 
derived parameters, ft,,,, h, and d\s. 

The key statistical observable discussed in this paper is the 
density fluctuation power spectrum P(k), the (ensemble-averaged) 
Fourier transform of the two-point density correlation function. 
In dimensionless form, the power spectrum can be written as 

k^Pik) 



2/;a — 



27r2 



(9) 



equivalent to the linear power spectrum in Eq. (O. Figure |7] 
shows the resulting dimensionless HaloFit power spectra for 
the 38 cosmological models at z = 1 (left panel) and at z = 
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Fig. 6. — Final simulation design: The five parameters under consideration 
are shown in red, projected onto two dimensions. The blue points show three 
derived parameters: f2,„, h, and rfj,. 
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Table 1 
Parameters for 38 Models 



# 










—w 


C^8 


# 






UJl, 


ris 


—w 










1296 


0.0224 


0.9700 


1.000 


0.8000 


19 





1279 


0.0232 


0.8629 


1.184 


0.6159 


1 





1539 


0.0231 


0.9468 


0.816 


0.8161 


20 





1290 


0.0220 


1.0242 


0.797 


0.7972 


2 





1460 


0.0227 


0.8952 


0.758 


0.8548 


21 





1335 


0.0221 


1.0371 


1.165 


0.6563 


3 





1324 


0.0235 


0.9984 


0.874 


0.8484 


22 





1505 


0.0225 


1.0500 


1.107 


0.7678 


4 





1381 


0.0227 


0.9339 


1.087 


0.7000 


23 





1211 


0.0220 


0.9016 


1.261 


0.6664 


5 





1358 


0.0216 


0.9726 


1.242 


0.8226 


24 





1302 


0.0226 


0.9532 


1.300 


0.6644 


6 





1516 


0.0229 


0.9145 


1.223 


0.6705 


25 





1494 


0.0217 


1.0113 


0.719 


0.7398 


7 





1268 


0.0223 


0.9210 


0.700 


0.7474 


26 





1347 


0.0232 


0.9081 


0.952 


0.7995 


8 





1448 


0.0223 


0.9855 


1.203 


0.8090 


27 





1369 


0.0224 


0.8500 


0.836 


0.7111 


9 





1392 


0.0234 


0.9790 


0.739 


0.6692 


28 





1527 


0.0222 


0.8694 


0.932 


0.8068 


1 



yj 




071 8 
u.uz i o 


U.O JUJ 


QQO 


7S^fi 
u. / J Ju 


7Q 
z.y 





1 7^fi 

iZJO 


0778 




Q1 ^ 


7087 


11 





1437 


0.0234 


0.8823 


1.126 


0.7276 


30 





1234 


0.0230 


0.8758 


0.777 


0.6739 


12 





1223 


0.0225 


1.0048 


0.971 


0.6271 


31 





1550 


0.0219 


0.9919 


1.068 


0.7041 


13 





1482 


0.0221 


0.9597 


0.855 


0.6508 


32 





1200 


0.0229 


0.9661 


1.048 


0.7556 


14 





1471 


0.0233 


1.0306 


1.010 


0.7075 


33 





1399 


0.0225 


1.0407 


1.147 


0.8645 


15 





1415 


0.0230 


1.0177 


1.281 


0.7692 


34 





1497 


0.0227 


0.9239 


1.000 


0.8734 


16 





1245 


0.0218 


0.9403 


1.145 


0.7437 


35 





1485 


0.0221 


0.9604 


0.853 


0.8822 


17 





1426 


0.0215 


0.9274 


0.893 


0.6865 


36 





1216 


0.0233 


0.9387 


0.706 


0.8911 


18 





1313 


0.0216 


0.8887 


1.029 


0.6440 


37 





1495 


0.0228 


1.0233 


1.294 


0.9000 



Note. — The five basic parameters for the 37 models that specify the simulation design and model which we use as an independent check. See text for definitions. 



(right panelfl Overall, the parameter space is well covered by 
these 38 models, the coverage being sufficient to accommodate 
upcoming weak lensing survey measurements. 




0.001 0.01 0.1 1 0.01 0.1 1 

k[Mpc"'] 

Fig. 7. — Dimensionless power spectra as given by HaloFit for the 38 
cosmologies specified in Table[T]at z = 1 (left panel) and z = (right panel). 



2.2. Emulation 

After specifying the design, the next task is to construct the 
emulator for predicting the nonlinear matter power spectrum 
within the parameter priors specified in the design. For an in- 
depth mathematical description of building such an emulator 
in the cos mological contex t we re fer the reader to iHabib et"al] 
dlooT) and Schneider et al.l 120081) . Here we explain the major 
ideas behind the process and show explicitly the emulator per- 
formance for our 37 model design (we exclude model due to 
the not very precise value for the Hubble parameter). As dis- 
cussed in the Introduction, we use HaloFit as a proxy for the 
full numerical simulations as a convenient foil to demonstrate 
and to test the overall procedure. 

'*Note that w e are using Mpc"' as units for k in this paper, different from 
IHeitmann et alj feOOS). We omit h since its value is different for each of the 38 
models we consider and it simplifies plotting the power spectrum. 



In order to construct the emulator, we model the simulation 
output using a /?,,-dimensional basis representation: 

r](k,z;9) = Y,^i(k,z)Wi(e) + e, 0e[O,lF^ (10) 

!=1 

Here, ri(k,z',0) represents the power spectra, over a 200 x 100 
grid of k and z values, with < z < 1 and -3 < iog^^k < 0.12. 
Over this support grid, the values ri{k,z',6) are determined by 
the five cosmological parameters denoted by 9. It turns out to 
be more convenient to model the power spectrum as 

than to work with A^{k,z) directly. This transformation re- 
duces the total dynamic range as well as enhances the baryon 
acoustic oscillation features in the power spectrum. The dimen- 
sionality refers to the number of orthogonal basis vectors 
{(piikjz), ■■■,<j>p^(k,z)}- As we will show later, Prj = 5 turns out 
to be an adequate choice for the present application. The pa- 
rameter pg is the dimensionality of our parameter space - with 
5 cosmological parameters we have pe = 5 (the two parame- 
ters being the same is a coincidence). As mentioned earlier, 
it is convenient to map the parameter ranges into [0,1] space. 
The Wi(d) are the weights of the basis vectors. The last term 
in Eq. (fTOl i. e, is the error term. Our main tasks in building the 
emulator are now: 

• Find a suitable set of orthogonal basis vectors (j)i{k, z). In 
our case, a principal components basis turns out to be an 
efficient representation, but alternative representations 
may be employed. 

• Model the weights w,(0). Here, we choose GP models; 
these have proven to be very well suited for representing 
functions that change smoothly as a function of 6, such 
as the power spectrum. 
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Table 2 

Derived Parameters for 38 Models 



# 






h 






# 






h 


dh 







0.2500 


0.7500 


0.7200 


14.24 


0.1 


19 


0.1940 


0.8060 


0.8120 


14.24 


<0.01 (7712) 


1 


0.4307 


0.5693 


0.5977 


13.59 


0.064 


20 


0.3109 


0.6891 


0.6442 


14.27 


0.15 


2 


0.4095 


0.5905 


0.5970 


13.80 


0.205 


21 


0.2312 


0.7688 


0.7601 


14.14 


<0.01 (21579) 


3 


0.2895 


0.7105 


0.6763 


14.10 


0.19 


22 


0.3317 


0.6683 


0.6736 


13.70 


<0.01 (11139) 


4 


0.2660 


0.7340 


0.7204 


13.99 


< 0.01 (5569) 


23 


0.1602 


0.8398 


0.8694 


14.48 


< 0.01 (4478) 


5 


0.2309 


0.7691 


0.7669 


14.11 


< 0.01 (2756) 


24 


0.1854 


0.8146 


0.8380 


14.21 


<0.01 (13138) 


6 


0.3059 


0.6941 


0.7040 


13.66 


<0.01 (19318) 


25 


0.4558 


0.5442 


0.5724 


13.76 


<0.01 (3033) 


7 


0.3310 


0.6690 


0.6189 


14.31 


0.225 


26 


0.2804 


0.7196 


0.6931 


14.05 


0.14 


8 


0.2780 


0.7220 


0.7218 


13.84 


< 0.01 (4320) 


27 


0.3357 


0.6643 


0.6387 


14.04 


0.08 


9 


0.3707 


0.6293 


0.6127 


13.93 


< 0.015 (2845) 


28 


0.3988 


0.6012 


0.6189 


13.66 


0.06 
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11 


0.2790 


0.7210 


0.7177 


13.82 


<0.01 (3928) 


30 


0.2810 


0.7190 


0.6626 


14.37 


0.155 


12 


0.2235 


0.7765 


0.7396 


14.43 


<0.01 (5901) 


31 


0.3791 


0.6209 


0.6394 


13.62 


<0.01 (17774) 


13 


0.3974 


0.6026 


0.6107 


13.77 


<0.01 (11549) 


32 


0.1922 


0.8078 


0.7901 


14.47 


0.115 


14 


0.3289 


0.6711 


0.6688 


13.74 


<0.01 (11803) 


33 


0.2634 


0.7366 


0.7286 


13.96 


<0.01 (2829) 


15 


0.2363 


0.7637 


0.7737 


13.89 


< 0.01 (9905) 


34 


0.3532 


0.6468 


0.6510 


13.71 


0.1 


16 


0.1981 


0.8019 


0.7929 


14.40 


0.025 


35 


0.3990 


0.6010 


0.6100 


13.77 


0.135 


17 


0.3586 


0.6414 


0.6305 


13.94 


<0.01 (5012) 


36 


0.2949 


0.7051 


0.6421 


14.41 


0.455 


18 


0.2578 


0.7422 


0.7136 


14.22 


<0.01 (5641) 


37 


0.2796 


0.7204 


0.7313 


13.71 


<0.01 (2971) 



Note. — The derived parameters, obtained from the basic parameters listed in TabIe[T]by applying the constraint on £/\. Only for model is the Hubble parameter 
picked by hand. The distance to last scattering is in Gpc, all other parameters are dimensionless. See text for details. 



2.2.1. The Principal Component Basis 

Before we determine the basis vectors to model the simula- 
tion output, we standardize the simulation power spectra in the 
following way. We first center the power spectra around their 
mean, given by ''Ij- resulting mean as a function of 

redshift z and wavenumber k is shown in the upper left corner 
of Figure [8] (Note that we divide A^(A;,z) by 2nk^^^ leading to 
the flattening of the power spectrum at high k.) After having 
centered the simulations around the mean, we scale the output 
by a single value so that its variance is 1 . 

The next step is the principal component analysis (PCA) to 
determine the orthogonal basis vectors (t>i{k,z) for modeling the 
simulation output following Eq. ( fTOl ). To carry out this step, 
we write the standardized power spectra as an n^, x m matrix, 
where n^; = 20000 is the number of support points for each 
power spectrum over the 200 x 100 k-z grid, and m = 37 is the 
number of simulations. The matrix then reads: 

. y. im. = r 77i;...;r?^7l. (12) 

Following iHabib et al.l ( 120071) . we apply a singular value de- 
composition to the simulation output matrix jsims giving 

y,^, = UDV\ (13) 
where U is an n^z x m orthogonal matrix, D is a diagonal m x m 
matrix holding the singular values, and V is an m x m orthonor- 
mal matrix. The PC basis matrix is now defined to be 
the first p,, columns of [UD/ ^Jm\ and the principal component 
weights are given by [^/mV]. 

In order to model the nonlinear matter power spectrum, 
we find that five principal components are sufficient to cap- 
ture all information. Therefore we have = 5 and = 
[0i;</'2;03;'?!'4;05]- The resulting five PC bases are shown in 
Figure |8] as a function of z and k. The fourth and the fifth com- 
ponents are already very flat - increasing the dimensionality 
further would not improve the quality of the data representa- 
tion. 



2.2.2. Gaussian Process Modeling 

The final step is to model the PC weight functions w,(6'), ; = 
1, . . . in Eq. (fTOl i conditioned on the known results from 
the 37 cosmological models. Gaussian process modeling is a 
nonparametric regression approach particularly well suited for 
interpolation of smooth functions over a parameter space. As 
mentioned previously, GP models are local interpolators and 
work well with space-filling sampling techniques. The GP (also 
called Gaussian random functions) is simply a generalization of 
the Gaussian probability distribution, extending the notion of a 
Gaussian distribution over scalar or vector random variables to 
function spaces. (For an excellent introduction to Gaussian pro- 
cesses, see Rasmussen & Williams 2006.) The Gaussian distri- 
bution is specified by a scalar mean /x or a mean vector and a 
covariance matrix - extending this to the GP leads to a specifi- 
cation of the GP by a mean function and a covariance function. 
We will take the mean function to be zero in the models de- 
scribed here. 

Before proceeding further, we briefly illuminate the concept 
of the GP with a simple one dimensional example shown in Fig- 
ure |9] The upper panel in Figure |9] shows three realizations for 
a mean-zero GP model (the data points have been normalized 
to a mean value of zero) of w{6) on 9\ ,...,9„ with n = 8: 



w(6i) = 



// \ 



( 



V w(f?8) ; V V / V 

or written more compactly: 
with 

Rij = ^x^{-\%-ejf}. 



(14) 

(15) 
(16) 
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mean PC basis 1 PC basis 2 




Fig. 8. — Mean (upper left corner) and five principal component (PC) bases derived from the output from the 38 HaloFit power spectra. The third axis shows 
the time evolution of the mean and the five principal components between redshift z = I and z = 0. While the first two principal component bases show significant 
variation, the fourth and fifth are already almost flat, indicating that the inclusion of higher principal components would not improve the quality of the emulator and 
the underlying GP model. 
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Fig. 9. — Example of a one-dimensional mean-zero GP model. The upper 
panels shows different draws from Eq. (TtJ. The middle and lower panel show 
an example of the conditioning of a GP model on two data points (blue dots). 
The data points could be measurements from an experiment or outputs from a 
simulation. The middle panel shows five realizations conditioned on the two 
data points. In this example, it is assumed that the data points are exact, so 
all realizations pass through them precisely. Our framework includes small 
uncertainties for the data points due to simulation noise and truncation eiTors 
from the finite basis set. Such uncertainties are denoted by A,, in the figure. 
The lower panel shows the mean of the curves from the middle panel. 



The symbol ^ means "distributed according to". The distance 
between locations 6*, and 6j is denoted by | |0, -6'j| |. The squared 
exponential form for the correlation ( fT6b is one among several 
possible choices. It enforces close to unity correlations between 
nearby w's. One can introduce a correlation length to control 
the range of the covariance, as we will do in the actual appli- 
cation below. The smoothness of the Gaussian process results 
from the infinitely differentiable nature of the correlation func- 
tion ( fTSI l. The overall covariance is determined by the multi- 
plicative prefactor A^' in Eqs. ( fT4b and ( fTSb . which in this ex- 
ample is set to unity. Restricting to these n = 8 locations, w has 



the probability density 

7r(w) = 27rT|E|"iexp|-iw^E"'w|, (17) 

where E = \~^R. The upper panel in Figure |9] shows three dif- 
ferent draws from this distribution. 

Suppose now that the value of the function w(6) is known 
(perhaps with some specified error) at a given set of 6** points. 
The aim then is to produce an ensemble of realizations of the 
Gaussian process consistent with the known data points at the 
specific values of 9* (passing through them if the error is zero), 
i.e., GP realizations conditioned on the data points. The second 
and third panel of Figure |9] show the example GP model condi- 
tioned on two data points (with zero error). The middle panel of 
Figure|9]shows five random realizations conditioned on the data 
points, while the lower panel displays the mean of the realiza- 
tions. The desired smooth interpolant w{9) is the mean over the 
ensemble, with the variance at any 6 providing a measure of the 
uncertainty in the interpolation. Note that as one approaches 
one of the specified 9* points this variance vanishes as 1^(6*) is 
specified exactly at those points. 

A naive approach to realize the conditioning would be to re- 
ject all the draws from Eq. ( [TtI i which do not go exactly through 
the data points. This is computationally very expensive and 
unnecessary. By specifying the data and the draws as jointly 
Gaussian random vectors [with separate means, but covariances 
assumed to be the same as specified in Eq. (fT4]i1. the conditional 
distribution becomes another known Gaussian process and may 
be sampled directly: 

w(02)|w(0r)~A^(S2iSnw(^r),S22-S2iS7;Ei2), (18) 

where 9i refers to the data points wher e w(9) is known exactly 
(see e.g. lRasmussen & Williamsll2006h . 

After this brief introductory interlude we now return to our 
full problem, where 9 lives in a p,, = 5 dimensional space and 
represents a cosmology with pe = 5 input parameters. We 
model each PC weight function Wi(9), i = 1,...,5 as a mean- 
zero GP 

Wi(9)^GP(Q,X;'R(^,0';p.,)). (19) 
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Here A,,,, is the marginal precision of the process and the corre- 
lation function is given by: 



Pe 



(20) 



This form is mathematically equivalent to that of ( fTSI ) - set 
P0 = 1, / = 1, and p = e~^l^ . The parameter p^^-^u controls the 
spatial range for the Zth input dimension of the process 
Under this parametrization, piv ,/ gives the correlation between 
WiiQ) and w,(6'') when the input conditions 6 and d' are iden- 
tical, except for a difference in 0.5 in the /th component. Our 
task is now to find X^, and pvv ,/ from the set of our simulations. 

From our 37 simulations, we first define a 5-component, 37- 
vector w, with ; = 1, ...,5: 

w*={wi{dX).-.^miy)\ (21) 

in comparison to Eq. ( fT4] i which is an 8-vector with only one 
component. The star indicates that we use our 37 input cos- 
mologies here and therefore the answer for is known at 
that point. Following Eq. ( fTsT i we assume that w* is normal- 
distributed with mean zero: 

w*~A^(0,S„,), (22) 
where S,, = diag(S,,, ...S,,^) and S,,,. = \}R(B*\p,,:). E„, is 
therefore controlled by five precision parameters A„ and the 
spatial correlation parameters held in p„. Next we have to spec- 
ify pri ors for each Ah., and for the p,, ;,/. Following iHabib et al.l 
(|2007|) . we choose r(flH ,^rt) distributions for the priors for A^v 
and beta(flp^, , bp^^, ) priors for the pn;// : 
1 



'^(A.v)=T;^^.v(^H.A„,)°"-'e 



7r(Pv,;,7) ■■ 



r(a,v) 
r(ap,, + fop„,) 



1,...,5, 



(23) 



p:;;r (1-P.v;-7A»-' /=l,...5j(24) 



r(flp„,)r(Z7p„,)^ 

with flu, = fo„, = 5, flp,,, = 1, and fep„, =0.1. The choices for fl^, 
and lead to a prior for A„ , of mean 1 and a prior standard 
deviation of 0.45 (The mean of a F distribution is given by a jh 
and the standard deviation by ^/a/W-). The choice of unit mean 
is consistent with the standardization of the GP for the w,-. 

The choices for flp^^ and bp^^. lead to a substantial prior mass 
near 1 [The mean of a beta distribution is given by a /(a+b) and 
the standard deviation by yJab/{{a + b)Ha + b+ 1)).] In general, 
the selection of these parameters depends on how many of the 
po inputs are expected to be active. 

Now we return to the actual information that we have, and 
from which we want to derive the weights w,-: the simulation 
outputs for the power spectra -q* for the 37 cosmologies. We 
arrange these outputs in an nt^m vector 

77* = vec([r,(0*);---;rK^3*7)])- (25) 
The simulation outputs have two sources of error: the error in- 
trinsic to the simulation (e.g., realization scatter, numerical er- 
ror) and the error due to the truncation in basis functions used 
to model rj{k\9) via Eq. ( fTOl i. We encapsulate the precision of 
the error in A,, (see middle panel of Figure |9] for an illustration 
of A,,) and we assume that the error e itself in Eq. ( fTOl i is in- 
dependent and identically normal distributed. We are now in a 
position to formulate the likelihood for 77*: 

L{i^*\w* , A„) cx A(7'^'/2 exp{-i A„(77* - ^w*f{Tj* - <i>w*)}, 

(26) 

where $ is a matrix composed from the 0,- basis vectors which 
we use to model the power spectra (see Eq. (fTOli). As for A„„ 



we specify the priors by a T{ar,,b^) distribution. We expect 
the data to be very informative about A,, and therefore choose 
the prior to be very broad with fl,, = 1 and = 0.0001. This 
prior allows for large values of A,, which force the GP model to 
nearly interpolate the simulation output. This will happen when 
the PC representation of the output is very good. 

This result is only an intermediate step, as our goal is to find 
the likelihood for the w, not for r; itself Fortunately, we can fac- 
torize Eq. ( |26] l to extract the likelihood for the weights easily. 
To do this we define vv as 



w = ($^$)-'$^?7*. 



(27) 



Note that w encapsulates the error due to the truncation of the 
base functions in modeling 77* . With this definition, it is easy to 
show that Eq. ( |26] l can be written as 

L(7?*|w*,A„) oc A('"^"»/2gxp{-iA^(w*-vv)^($^$)(w*-vv)} 



y \ (m(ni-pn)/2 

xexp{-^A^77*^(/-$($^$)-'$^)7f}, 



(28) 



with m = 37, = 5, and iii^- denoting the k-z support for each 
power spectrum is measured. Note that in the first line of 
Eq. ( l28T l, w is completely separated from the rest of the like- 
lihood expression. We can use this factorization to represent 
the likelihood in a dimension-reduced form: 

L(vv|w*, A„) cx A('"""»/2exp{-iA„(w*-w)^($^$)(w*-vv)}, 

(29) 

where the remaining terms from Eq. ( |28] ) are absorbed in a re- 
defined Gamma distribution prior for A^, T{a[^,b'^) with 

m(nkz-pri) 



a^ = a^ + - 2 
1 



b'^ = br, + -v^(I-H^^<^>r^^)V. 



(30) 



(31) 



It is useful to recap what has been done so far: We began 
with the normal likelihood for 77 with the Gamma distribution 
prior for A,, : 

77*|w*,A^-A?($w*,A-^/„,J, A^-r(fl,„/7^). (32) 



o 



o 



o 



Fig. 10. — Boxplots of posterior samples for each p^-ji for the nonlinear 
matter power spectrum. Boxplots offer a convenient way of showing the dis- 
tribution of data using just five numbers. The blue box itself contains 50% of 
the data, the lower edge indicating the 25th percentile and the upper edge, the 
75th percentile. The red (center) line denotes the median. If the red line is 
not at the center of the box, the data is skewed. The black lines (or whiskers) 
extend out to the full range of the data. With our parametrization, a box value 
close to 1 indicates that the parameter is inactive, i.e., the PC is not changing 
much under the variation of that parameter 
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Fig . 11. — Posterior mean estimates of the principal component weight functions wj {8) to n'5 (8). Here the prediction points 6 are over a grid of ui,,, and «s values, 
while the remaining parameters are held fixed at their midpoints. The cosmological parameters are displayed in the normalized [0,1] space. 



Due to the relation 

L(r]*\w*,Xr,) X 7r(Xrj;an,br,) oc L(vv|m',A^) (33) 

we can derive the likelihood for w: 

vv|w*,A^^A^(w*,(A„$^$)-'), X^^Tia'^X,). (34) 
Next, w* is integrated out, leading to the posterior distribution 
|vv) cx (35) 



|(A„$^$)-' + E„ 



X exp{--vv^([A„$^$]-' + E,,)"'vv} 



(36) 



X Xrl^ 



a' -I 



Pr, 



Pri Pe 



ho-\ 



As detailed in iHabib et alj (l2007h this posterior distribution is a 
milepost on the way to creating the emulator for the power spec- 
trum. It can be explored via MCMC and contains much use- 
ful information about the parametric dependence of the power 
spectrum, as derived from the numerical simulation results at 
the finite number of design points. 

For example, results for the p„ are shown in Figure [10] in 
form of boxplots. With our definitions (see also Habib et aL. 
l2007h . an input / is inactive for PC / if pw,ii = 1- In the one 
dimensional example, shown in Figure |9] the level of activity 
can be roughly estimated by how strongly the curve varies. An 
almost flat curve would lead to /?„, ~ 1 . Figure [10] shows that 
ujm and (78 are clearly the most active parameters influencing the 
power spectrum. The equation of state parameter w is also very 
active. This is mostly due to the fact that the Hubble parameter 
is changing for each cosmology and the best-fit value for h for 
each model is used. This in turn influences the behavior of the 
power spectrum with varying equation of state parameter w. 

The last step for building the complete emulator is to draw 
from the posterior distribution ( [35l l at any given 6. We consider 
the joint distribution of w and a predicted weight We at a new 



input parameter setting 6^: 



N 



(A„$^$)-i 











+ S,v.h,,(Ak', Pw) 



(37) 

where E„,.,j,^ is obtained by applying the prior covariance rule 
to the augmented input settings that include the original design 
and the new input setting (6e)- Following Eq. ( [TSl l we have 

(38) 



We W ' 



'NiV2iV[lw,V22-V2iV[IVn), 



where 



V = 



Vn V12 
V21 V22 







(39) 



is a function of the parameters produced by the MCMC output. 
Hence for each posterior realization of A,, , A„, , p^, a realization 
of w can be produced and the emulator is completed. 

Figure [TT] shows the results for the PC weights wi to W5 cor- 
responding to the PCs shown in Figure [8] We show the results 
as a function of two of the active parameters, u!„ and «j, while 
holding the remaining 3 parameters fixed at their midpoints. 
Note the very smooth behavior of the weights as a function of 
the parameters. 

2.2.3. Emulator Performance 

In order to evaluate the accuracy of the emulator we generate 
a second set of power spectra with HaloFit within the prior 
parameter ranges. For this set we choose the input cosmologies 
randomly, but still insuring that the constraint on the Hubble 
parameter is satisfied. We then predict the results for those cos- 
mologies with the emulator scheme and compare them to the 
HaloFit output, the "truth" in this case. The results for the 
residuals are shown in Figure [T2]for three redshifts, z = 0, 0.5, 
and 1 . The middle 50% of the residuals (dark gray band) are 
accurate to 0.5% or better over the entire A;-range and for all 
three redshifts. All predictions have errors less than 1%. This 
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result shows that a simulation set with as small a number as 37 
cosmologies is sufficient to produce a power spectrum emulator 

accurate at 1%. 

In Habib et al.l (l200l several other convergence tests are dis- 
cussed and demonstrated. These tests show that emulator per- 
formance can improve considerably - by an order of magnitude 
- if either the number of simulation training runs is increased 
or the parameter space under consideration is narrowed. In the 
present paper we follow the second strategy, restricting the pri- 
ors as much as is sensible given the current and near-future ob- 
servational situation. 



z=1 z = .5 z = 
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Fig. 12. — Emulator performance at three redshifts, z = I, 0.5 and (left 
to right). The emulator is tested on 10 additional HaloFit runs within the 
parameter priors. The emulator eiTor with respect to the HaloFit results is 
shown. The central gray region contains the middle 50% of the residuals, the 
wider light gray region, the middle 90%. Errors are at the sub-percent level. 



2.3. Sensitivity Analysis 

After the emulator has been built it can be used to explore the 
behavior of the power spectrum within the parameter priors in 
more detail. An obvious use of emulation is to carry out a sen- 
sitivity analysis, i.e., study the behavior in the power spectrum 
as the underlying cosmological input parameters are varied. 




-3 -2 -1 -3 -2 -1 -3 -2 -1 

log k log k log k 



Fig . 13. — Sensitivity analysis for each of the five parameters at redshift z = 
0, 0.5, and 1 (bottom to top). The 3'-axis shows the deviation of the log of the 
power spectrum from the nominal spectrum where each parameter is set at its 
midpoint. The light to dark lines correspond to the smallest parameter setting 
to the largest. Due to the tight constraints on uji, from CMB measurements, 
which led us to choose a rather narrow prior, ui/, variation leads to very little 
change in the power spectrum. 



An example of this is represented in Figure [T3] where varia- 
tions of the power spectrum at three redshifts z = Q, z = 0.5, and 
z = 1 are shown. The reference power spectrum is the power 
spectrum obtained if every parameter is fixed at the midpoint of 
its prior range, i.e., in this case, for the cosmology uj,„ = 0.1375, 
ojh = 0.02215, n.5 = 0.95, w = -1, = 0.758. (This power spec- 
trum is close to the mean of the 37 models from our design but 
not the same.) Next, only one parameter is varied between its 
maximum and minimum value while the other four parameters 
are fixed at their midpoints. In Figure [T3] from left to right we 
vary oj,,,, uib, n^, a^, and w, showing the difference between nat- 
ural logarithm of these two power spectra. We note that the 
Hubble parameter is different for each power spectrum shown 
in the figure since it is separately optimized for each cosmology. 

The results contain information about the scales at which the 
power spectrum is most sensitive to each parameter and about 
parameter degeneracies. For example, it is clear that the power 
spectrum is relatively insensitive to uji, at any scale or redshift 
and therefore will be difficult to constrain from power spectrum 
measurements alone. In the quasi-linear to nonlinear regime 
at A: ^ 0.1 - 1 /iMpc"', the power spectrum holds significant 
information regarding erg and w, but degeneracies become an 
issue. Very large scales are particularly sensitive to the spectral 
index and uo„,. 

The sensitivity analysis is also helpful in the targeted aug- 
mentation of simulation designs. If the accuracy of the emula- 
tor is not sufficient for the problem of interest, one would like 
to improve it by adding additional simulations. These simula- 
tions would then involve variations of the most active parame- 
ters while keeping the other parameters more or less fixed. The 
augmentation of existing designs is an active field of research 
in statistics, and potentially important in precision cosmology 
applications. 

3. CONCLUSIONS 

The last three decades have witnessed unprecedented progress 
in cosmology. From order of magnitude and factor of two es- 
timates for cosmological parameters, we now have measure- 
ments at 10% accuracy or better These measurements have 
revealed one of the biggest mysteries in physics today: a dark 
energy leading to the acceleration of the expansion of the Uni- 
verse. In order to understand the origin, nature, and dynamics 
of this dark energy - or to prove that the acceleration is due to a 
modification of gravity on the largest length scales - the accu- 
racy of the measurements must be further improved. The next 
step, as defined by near-term and next-generation surveys is to 
obtain measurements at the 1 % accuracy level. This puts con- 
siderable stress on the quality of theoretical predictions, which 
have to be at least as accurate. Three major probes of dark 
energy - baryon acoustic oscillations, weak lensing, and clus- 
ters - are based on measurements of the large scale structure 
in the Universe. In order to obtain precise predictions for these 
probes, expensive, nonlinear simulations have to be carried out 
and ways must be found to extract the needed information from 
a limited number of such simulations. 

In this paper, we have demonstrated that if very accurate sim- 
ulations are available, 1 % accurate prediction schemes can be 
produced from just tens of high-accuracy simulations. The fo- 
cus of this paper is the nonlinear matter fluctuation power spec- 
trum, but the general scheme applies to any other cosmological 
statistic, e.g., the halo mass function. 

We have introduced a set of 38 cosmological simulations, the 
Coyote Universe suite, all of which satisfy the 1 % error con- 
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trol criterion dHeitmann et alj|2008l) for the power spectrum up 
to A; ^ 1 ZiMpc"'. Following the results of the current paper 
(demonstrated on HaloFit), from these simulations we can 
generate an emulator for the nonlinear power spectrum which 
has essentially the same accuracy as the simulations them- 
selves. The high accuracy attained from a small number of 
simulation inputs is due to (i) an interpolation method based on 
a sophisticated simulation design and GP modeling which has 
been developed and refined in the statistics community over the 
last decade to address problems of the nature described here, 
and (ii) the excellent parameter constraints from CMB measure- 
ments, which allow us to base our emulator on relatively narrow 
parameter priors and therefore ease the interpolation task. 

This paper is the second in a series of three papers with the 
final goal to provide a high-precision emulator for the nonlin- 
ear power spectrum out to ^ 1 /iMpc"'. The first paper of the 
series (Heitmann et al. 2008) demonstrated that the required ac- 
curacy to predict the power spectrum from A^-body simulations 
can be achieved. The current paper introduces the cosmolo- 
gies underlying the Coyote Universe simulation suite, explain- 
ing and demonstrating success of the emulation technology us- 
ing HaloFit as a proxy for the simulation results. Our pre- 
diction scheme can achieve 1% accuracy from only a limited 
number of simulations: approximately 37 cosmological models 
are adequate for this purpose. The third and final paper will 
present results from the simulation suite discussed in this paper 
and will include a power spectrum emulator built around them. 
This emulator will be publicly released. 
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supported in part by NASA and the DOE. We would like to 
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